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We present results from direct numerical simulations of the 
Kraichnan model for passive scalar advection by a rapidly- 
varying random scaling velocity field for intermediate values 
of the velocity scaling exponent. These results are compared 
with the scaling exponents predicted for this model by Kraich- 
nan. Further, we test the recently proposed fusion rules which 
govern the scaling properties of multi-point correlations, and 
present results on the linearity of the conditional statistics of 
the Laplacian operator on the scalar field. 
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In the last equation the scaling of h{R) is expressed nor- 
malized with respect to C, the outer scale of the velocity 
field. For physically realisable fields may vary between 
and 2. 

Our aim is to express the statistical properties of the 
scalar field in terms of the parameter (^h- The statistics 
is characterized by the n-point correlators, defined as 
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As one of the simplest realisations of a model with 
turbulent statistics with non-trivial scaling exponents, 
the Kraichnan model of advection by a white-in-time 
scaling velocity field has attracted much recent attention 
|^|-^. The model is analytically tractable, in the sense 
that its statistical description may be reduced to a set 
of closed form differential equations for the n-order cor- 
relation functions. The model concerns the equation of 
motion for a passively-advected scalar field T driven by 
a velocity field u: 



d_ 
dt 



T{x, t) + u{x, t) ■ VT{x, t) = kV'^T{x, t) -f f{x, t), 



(1) 



where k is the molecular diffusivity. The velocity field 
is taken to be a Gaussian, white-in-time, incompressible 
homogeneous scaling random field. Statistical stationar- 
ity is achieved through the forcing /, which is also taken 
to be delta-correlated in time, statistically homogeneous 
and isotropic, and to exhibit only large-scale spatial com- 
ponents. The parameter of interest in this model is the 
scaling exponent characterizing the so-called eddy dif- 
fusivity tensor hij(R) which contains the relevant infor- 
mation about the random velocity field u{r, t): 

hij{R)= J dT(^[u,{r + R,t + T) -u,{r,t + T)] 

x[uj{r + R,t)^Uj{r,t)]) . (2) 



The notation (• • •) refers to ensemble averaging. Under 
the conditions that the velocity field exhibits fast tem- 
poral decorrelation, scaling and incompressibility, hij{R) 



takes the d-dimensional form 0] 



hij{R) = h{R) 



Ch + d-l 
d-1 ' 



C/i RiRj 
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One expects the correlators to be homogeneous func- 
tions of their arguments, J>i(Ari, Ar2, .., Ar„) ~ 
A^"^„(ri, r2, .., T„), and one hopes to determine the de- 
pendence of the scaling exponents Cn on (^h- In this 
model, the rapid temporal decorrelation of the velocity 
allows one to derive a set of closed equations for these 
correlation functions [|] 
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where J-2n is a function of the 2n variables ri, r2, r2„ 
and J-2n~2 is a function of the 2n — 2 variables 
ri, r2, r2n except for and rp. <i>o is the forcing 
correlation and may be eliminated using the two-point 
equation. Only the 2rtth order moments are considered 
as by isotropy odd moments vanish. For n = 1 these 
equations are readily solvable, leading to the exact result 



C2 = 2 - 0. . 



(7) 



For n > 2 the equations are difficult to solve analyti- 
cally for arbitrary values of C^, and to date only certain 
limits have been treated. The limit of k — > 0, C/i ^ 
(in that order) has been examined perturbatively in [Q. 
This limit is not realisable in direct numerical simulations 
due to numerical instabilities caused by small diffusivi- 
ties; moreover fields with scaling exponents approaching 
zero become increasingly spatially rough and are very 
difficult to produce and treat reliably numerically. In 
Q the perturbative small parameter was C,h/d, with d 
the spatial dimension; which requires either the difficult 
Qh ^ Q limit or the numerically inaccessible case of large 
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dimension. The regime of C/i ^ 2 has also been treated 
perturbatively in The only theory which treats the 
intermediate span of physical fields requires a closure that 
is not rigorous [^J^ , and it is with the prediction arising 
from this theory that we will be able to make a compar- 
ison. Further we test in detail some of the more general 
scaling predictions afforded by the fusion rules for fluid 
dynamics developed in and the particular statistical 
assumptions with respect to conditional statistics utilised 
in the theory of in obtaining predictions for the scal- 
ing exponents. 

The crucial assumption arises in the context of the 
equation for the nth order structure functions, defined 
Sn{R) = {{Tix + R)-Tix)r): 

R'~''^R'''HR)4B^2niR) = J2n{R) ■ (8) 

oR oR 

The function J2n{R) derives from the dissipative term 
and is given by 



(9) 



where SfjT{x) = T{x + R) — T{x). One may determine 
directly that J2{R) = 4e, the mean dissipation (indepen- 
dent of R). 

In order to obtain the scaling exponents Cn of the nth 
order structure functions, one needs to evaluate J2n{R)- 
In light of (^ and the exact result (0) one sees that J2„ 
must have a scaling form which agrees with 

J2n{R) = nC2nJ2S2n{R)/S2{R) for n > 1 . (10) 

This result can be derived without reference to (|^) using 
the fusion rules derived in In either way the coeffi- 
cients C2n are undetermined. Kraichnan proposed that 
C2n = 1 for all n. In this case one obtains from (^) a 
quadratic equation determining the CnS: 
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C2-d+ V(C2+d)2+4C2(n-l) 



(11) 



As has been pointed out in Q this assumption bears a 
strong relation to the conditional statistics of the Lapla- 
cian of the field. One may rewrite J2n{R) in terms of the 
average of the Laplacian conditioned on the value of a 
difference of T across the length scale i?, 5bT{x): 



J2n{R) = -2nK J d5RTP{dRT)[SRT]' 
X {V^T{x)\SrT{x)) , 



(12) 



One way to ensure that J2n {R) has the scaling ( |To| ) is for 
the conditional average to satisfy 

{\/^T{x)\6rT{x)) = Ce6RT{x)/KS2{R). (13) 

Hence a linear behaviour of the conditional average of the 
Laplacian is intimately connected with the determination 
of the scaling exponents. 



The model has been studied by direct numerical sim- 
ulations in 1^ with Qi = 1. These simulations have been 
criticised for the method of generation of the velocity 
field; two fixed scaling fields were swept past each other 
in orthogonal directions at a constant rate. In doing so 
one may lose isotropy in a way that can influence the ap- 
parent numerical values of the measured exponents. In 
our simulations we have evolved a scalar field in two di- 
mensions on a 1024^ grid. The scaling velocity field was 
implemented by Fourier transforming a set of fc-vector co- 
efficients which were each chosen randomly from a Gaus- 
sian distribution scaled to a standard deviation propor- 
tional to The direction of the fcth component 
Uk was chosen such that k ■ Uk = 0- To reduce computa- 
tion we have used an isotropised version of the method 
employed in namely we generate two fixed realisa- 
tions and shift them with respect to one another in order 
to obtain rapid variation. At each time step the two 
fields are independently shifted by a step of random size 
and direction. The fields are renewed after around every 
500 time steps to reduce any temporal correlation that 
this method might induce. We checked that the results 
are insensitive to a more frequent refreshment of these 
fields. The spatial discretisation is second order, and 
the time evolution was performed using an explicit Eu- 
ler scheme. The forcing was implemented by stimulating 
at every time step one of the nine smallest wavenumbers 
with an amplitude chosen from a Gaussian distribution. 
Our initial conditions for the scalar field (for a given value 
of C/i) were Gaussian random with the 2nd order scaling 
exponent distinct from the expected result of 2 — C^, and 
truncated in k space. Typically, saturation to statisti- 
cal steady state required about thirty million time steps 
on the CRAY J90. We have converged results for three 
values of (h, i-e. 0.6, 1.0 and 1.2. The diffusivity in ev- 
ery run was chosen to obtain the longest possible inertial 
range while retaining stability in the small scales. 

In Fig. |l| we present a typical realization of the scalar 
field for = 1.0. It shows significant development of 
small scale structures. In Fig. ^ wc present the structure 
functions S'„(i?) as a function of R for the three values 
of Ch, computed using spatial averaging over single re- 
alizations after statistical stationarity was reached, and 
then time averaging over one hundred snapshots taken at 
intervals of ten thousand time steps. This figure shows 
that we have one and a half decades of scaling, or "iner- 
tial range" . 

Figs. H displays the dependence of Cn on n for the three 
values of (^h- Also shown is the prediction of Kraichnan 
for these values. It is evident that for the three parame- 
ter values tested we have close agreement. In the figures 
we display also the odd values for the exponents. These 
were calculated from the field by taking absolute values; 
strictly this is not covered by the theory but one sees 
here that they smoothly interpolate the law for the even 
orders. We remark that although the grid is relatively 
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small the structure functions display well-developed scal- 
ing ranges for orders as high as 12. The relatively good 
statistics resulted from averaging over many snapshots. 
We checked however that also the single-time realisations 
appear to be well self-averaged. 

Note that for (^h — 1.0 the agreement between the nu- 
merically computed value of ^2 and Eq.(|^) is best. We 
believe that the reason for this is simply due to the diffi- 
culty of creating a velocity field with precise scaling on a 
finite grid. It is interesting that in fact the scaling in the 
passive scalar field appears cleaner than that which can 
be obtained by the Fourier transform method described 
above in grids of this size. If we check our apparent real 
space scaling exponents for the velocity field we find that 
the minimum error between the input (h in A:-space and 
the observed one occurs precisely at C/i = 1- However 
the higher order scaling exponents do not seem to be as 
sensitive to this discrepancy. 

The quality of the prediction (|ll] ) can be independently 
tested by verifying that the coefficients C„ are close to 
unity, and that the conditional average (^3|) is indeed 
linear with the right i?-dependent prefactor. To this end 
we computed from the simulation the quantities J„ (i?) of 
Eq.(^. J2 was confirmed to be constant throughout the 
inertial range. In Fig. ||we present J„(i?) as a function of 
nJ2Sn{R)/2S2{R) for n = 2,4, 6, 8, 10 and inertial range 
R. The dashed line is the line y — x, and we see that it 
passes through the data without any adjusted parameter. 
The coefficients C„ were obtained from the data for a 
range of values of R and n, and were found to be very 
close to unity in the inertial range, see inset in Fig. ^. 

Finally we can check the postulated linearity of the 
conditional average (13). These quantities were calcu- 



lated for a range of R values in the inertial range by 
averaging over several directions of R. The results are 
displayed in Fig. [Hj 

Our conclusions from these simulations are that the 
postulates that lead to the prediction (^Tj) for the scaling 
exponents (i.e linear conditional averages, C2n = 1) are 
very well supported by the numerical data. As a result it 
is no surprise that the measured scaling exponents agree 
very closely with their predicted values. Due to the lim- 
itations of the computational techniques one cannot of 
course state that precise agreement is observed. It is our 
conviction however that the conditional average is very 
close to being linear; a persistent failure to prove the lin- 
earity mathematically may indicate that this property is 
not exact. It seems however very worthwhile to probe 
this question further to understand the close agreement 
between simulations and (O). 
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FIG. 1. Typical realization of the scalar field 



FIG. 2. Log- log plot of the structure functions S'n(-R) as a 
function of R for n = 2, 4, 6, 8, 10. 



FIG. 3. The scaling exponents as a function of n = 2 — 10 
for three values of (h- The numerical data (error bars) are 
compared to the analytic prediction Eq.(ll) (dotted line) 



FIG. 4. Jn{R) as a function of the fusion rule prediction 
nJ2S„{R)/2S2{R) with C„ = 1 for Ch = 1-2. An independent 
measurement of C„ is exhibited in the inset. The other values 
of show equivalently good agreement. 



FIG. 5. Conditional averages normalised by the scaling of 
Eq.(13) calculated for the field with (^h ~ 1-0, and from a 
single realization. Equally satisfactory results were obtained 
for the two other values of ^h- 
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